/***************************
Make the interpolation in 2D (sort of..)


*****************************/

using namespace std;

extern vector<double> tab_fuv;
//time is the x-axis
//mass is the y-axis
double interpolate_imf_fuv( double time_ind, double mass_ind, 
		     long size_x) {

  

  //initialize variables
  long fx, cx, fy, cy; 
  double fracx, fracy;
  double p1, p2, p3, p4;
  double interp_mass1, interp_mass2;
  double pval;
   
  
  //definte some frequently used variables
  fx=(long) floor(time_ind);
  cx=(long) ceil(time_ind);
  fracx=time_ind-fx;
  fy=(long) floor(mass_ind);
  cy=(long) ceil(mass_ind);
  fracy=mass_ind-fy;
 

  long sfy=size_x*fy;
  long scy=size_x*cy;
  
 
  //find the values of the 4 points around our
  //value... specifically p1=upper left, p2=upper right, 
  //p3=lower right, p4=lower left
  p1=tab_fuv[fx+scy];
  p2=tab_fuv[cx+scy];
  p3=tab_fuv[cx+sfy];
  p4=tab_fuv[fx+sfy];
//	cout<<"--"<<fx<<"\t"<<cx<<"\t"<<scy<<"\t"<<sfy<<endl;
// 	cout<<p1<<"\t"<<p2<<"\t"<<p3<<"\t"<<p4<<endl; 
//	cout<<"fracy fracx  "<<fracy<<"\t"<<fracx<<endl;
  //interpolate along the mass direction
  //This is only a local thing that doesn't know about the convention of 2-1 (RdS)
  interp_mass1=(p1-p4)*fracy+p4;
  interp_mass2=(p2-p3)*fracy+p3;
//	cout<<"interp_masses  "<<interp_mass1<<"\t"<<interp_mass2<<endl;  
  //now interpolate along the time direction
  pval=(interp_mass2-interp_mass1)*fracx+interp_mass1;
  
// cout<<pval<<endl;
  //now we return the 10^log(variable)
  pval=pow(10, pval);
  return pval;
}
